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Different phases in open driven systems are governed by either shocks or rarefaction waves. 
A presence of an isolated umbilic point in bidirectional systems of interacting particles sta- 
bilizes an unusual large scale excitation, an umbilic shock (U-shock) . We show that in open 
systems the U-shock governs a large portion of phase space, and drives a new discontinuous 
transition between the two rarefaction-controlled phases. This is in contrast with strictly 
hyperbolic case where such a transition is always continuous. Also, we describe another 
robust phase which takes place of the phase governed by the U-shock, if the umbilic point is 
not isolated. 

I. INTRODUCTION 

Many intrinsically nonequilibrium phenomena can be observed already in simplest systems of 



driven diffusing particles 



2|] , which are paradigmatic models of systems far from equilibrium and 



find a wide range of applications in biological, social and physical contexts [3-jJ. Driving forces 
due to bulk fields or boundary gradients lead to steady state currents that invalidate the condition 
of detailed balance and give rise to remarkable features which have no equilibrium counterparts, 
such as boundary driven phase transitions, spontaneous symmetry breaking and hysteresis in one 
spatial dimension. Models with two or more conserved species of particles exhibit particularly rich 



fl. 



behaviour 

The evolution of driven systems on large spatio-temporal scales is governed by two fundamental 
types of excitations: shocks, which carry discontinuities, and rarefaction waves, which are con- 
tinuous self-similar solutions of the hydrodynamic limit equations (3). Various properties of the 
fundamental excitations like stability, speed and morphology are determined by a scalar or vector 
function which relates steady macroscopic currents to average particle densities, the so-called cur- 
rent density relation. The topology of the current-density function such as the number of extrema 
and saddle points determines qualitative features of the large scale dynamics and in particular 
the number and character of the different stationary phases and phase transitions that one can 



observe in the underlying microscopic model 



8l4l0l]. In this way microscopic details of local par- 
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tide interactions are largely irrelevant as long as they produce a certain type of a current density 
relation. 

It was noted ll| that bidirectional particle systems, in which bulk hopping rates of oppositely 



moving interacting species possess the left-right symmetry, have a special property: their current- 
density function has an umbilic point. A generic umbilic point is a point on a current-density 
surface where the two characteristic velocities coincide, which breaks the usual assumption of 
strict hyperbolicity 12j. For bidirectional systems with left-right symmetry, both characteristic 
velocities vanish at the umbilic point which can be isolated or not, depending on a strength of 
an interaction between the species. An isolated umbilic point, in an open system under maximal 
feeding regime (a regime where particles enter and exit the system freely) stabilize a large scale 
excitation reminiscent of a shock wave, but which should be unstable according to usual shock 
s t ahil ity criteria Q,Q. The new excitation was caned an n mbill e shoe,, or a U-shoeh, and 
studied on microscopic scale fll]. 

In present article we determine a domain on the phase diagram which is controlled by a U- 
shock, and identify a boundary driven phase transition that it governs. Boundary-driven phase 
transitions in driven systems, which are caused by adiabatical changes of boundary conditions have 
no equilibrium counterparts 15], [s], [if]]. They may be continuous or discontinuous, depending on 
whether the order parameter changes across the transition in a continuous or discontinuous way. We 
show that a U-shock governs a discontinuous phase transition from one rarefaction-wave controlled 
state to another. Such a transition in a usually considered strictly hyperbolic systems (without 
an umbilic point) is always continuous If the umbilic point is not isolated, the U-shock is 

no longer stable, but on its place we find another robust phase, a homogeneous bulk state with 
densities matching the umbilic point. This state, which we call the umbilic state, has the same 
stability domain as the U-shock phase. A boundary driven phase transition leading to the umbilic 
state occurs via a continuous transition. We expect our resuts to be generally valid for any system 
with an umbilic point; however for definiteness we consider a model for which the current-density 
relation is known exactly. 

The plan of the paper is the following: In Sec. [II] we introduce our model. In Sec. IHII wc 
discuss splitting of the physical region according to signs of characteristic velocities, and review 
the U-shock and the umbilic phase. In Sec JI VI we describe phase transitions from and to the phases 
controlled by umbilic point, along a trajectory, where boundary rates are changed adiabatically. 
There, we also describe the domain of stability of the umbilic point- controlled phases. We finish 
with conclusions and perspectives. Appendices contain necessary technical details. 
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Figure 1: Bidirectional two-chain model. For solvability, the rates must satisfy a = /? = e = l,7 = e l/ where 
v is the interchain interaction constant 19J. Coupling to boundary reservoirs is indicated by boxes marked 



L (the left reservoir) and R (the right reservoir). 

II. THE BIDIRECTIONAL MODEL WITH BOUNDARY RESERVOIRS 

Our model describes particles with repulsive hard-core interaction which hop unidirectionally 
along two chains of N sites: One chain for right-hopping particles and another chain for left- 
hopping particles. At each instant of time the system is fully described by occupation numbers 
rik £ {0, 1} (for the right movers) and £ {0, 1} (for the left-movers). A right-moving particle 
at site k can hop to its neighbouring site k + 1 provided it is empty, with a rate that depends on 
the occupancies at sites k, k + 1 on the adjacent chain, see Fig{TJ E.g. a particle hops with rate (3 
if the adjacent sites are both occupied, etc. For clarity of presentation and analytic simplification 
we shall keep only one rate 7 = e u different from others, setting all remaining rates to 1, 



a = fi 



7 



(1) 



Then the parameter 



= 7-1 



(2) 



which ranges from —1 to 00, measures the interaction strength between the left- and right-moving 
species. For Q = the model reduces to two independently running totally asymmetric exclusion 
processes 

The bulk dynamics of particles (see FigQ]) is complemented with boundary conditions: we 
consider open boundaries where at the left end of the chain a right mover can enter the chain and 
it can leave it at the right end. Left movers are hopping to the left with the same dynamic rules. 
The boundary hopping rates are chosen so as to correspond to particle reservoirs with effective 
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densities of right- and left movers ul,vl at the left boundary, and ur,vr at the right boundary, see 
Appendix A. Note that generically ^ vl,ur ^ vr , so that while the bulk dynamics is left-right 
symmetric, the entrance and exit rates for different species are not. After a certain transition 
period, the system will approach a stationary state, characteristics of which (the average flux, the 
density profile, the correlations) do not depend on time. 

For our dynamical Monte-Carlo simulation we choose the following random sequential update 
procedure. For a chain of length N, i.e. a system of 2N sites (numbered i = 1,2, ...N for right 
movers and i = N + 1, N + 2,...2N for left movers) one Monte-Carlo step consists of 2N + 2 
uniform drawings of an integer random number w in the range < w < 2N + 1. If < w < N, the 
configuration of right movers is updated. If w = 0, and the left boundary site i = 1 is empty, we 
fill it with a particle with a rate parametrized by boundary reservoir densities ul,vl, see Appendix 
El If w = N and the respective site contains a particle, we remove it with a rate parametrized 
by boundary reservoir densities ur, vr, see Appendix [A] For intermediate < w < N, if site w 
contains a particle, a hopping is performed on the right neighbouring site with given rates (pQ), 
provided it was empty. The update of the left movers is done analogously. We start from an empty 
lattice and after a transient time we measure site occupancies , , and take averages over many 
Monte Carlo steps and many histories. We choose a system size up to L = 500 sites in each chain. 
The transient time for L = 500 is 10 5 Monte Carlo steps, and averaging over up to 10 histories is 
done. 

In contrast to our study U| which was focused on microscopic features, here we focus on large- 
scale hydrodynamic behaviour of an open system with an umbilic point. To this end, we also use 
an improved meanfield approach, described in Appendix |Bj Results obtained by the stochastic 
approach and the meanfield approach agree well, both for system dynamics and for steady state 
global averages, due to product- measure steady state property on a ring (|A1|) . 



III. SPLITTING OF PHYSICAL REGIONS ACCORDING TO SIGNS OF 
CHARACTERISTIC VELOCITIES 

Characteristic velocities c\(u, v) and 02(11, v) are velocities with which infinitesimal perturbations 
are propagating, on top of a stationary homogeneous background with average densities of right- 
and left-moving particles u and v. As such, they play a fundamental role in the stability of large 



scale excitations 



ISO- 



The characteristic velocities can be obtained by solving an eigenvalue problem for a flux Jacobian 
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(Dj)* k = c fc ^ fc where 



dji 8h 

IS Si 

du dv 



For our model (pQ) the particle currents ji and ji can be obtained analytically, see [lij for details, 
and are given by 

jt(u, v) = u(l -u) + QVtn{u, v)Q 00 (u, v) (4) 
h(u,v) = -ji(v,u) = -v(l -v)- Qtlii(v,u)Q 00 (v,u), (5) 

where Qn and S7oo are are stationary probabilities to have two adjacent particles and two adjacent 
holes, 



(u + v-l)Q-l + ^((u + v-l)Q- l) 2 + AQuv 
nu = 2Q (6) 

r^oo = i — u — v — fin. 

Due to the hardcore exclusion, the average densities of the right and left-moving particles may 
only take values between and 1. The whole physical region of < u, v < 1 is then splitted into 
regions G CTT with different signs of characteristic velocities, which is illustrated in Figl2j Subscripts 
a, t denote signs the characteristic speeds c\ and C2, i.e. a = 0,+,— correspond to c\ = 0,ci > 

0, c\ < 0, and similarly for r. E.g. we name by G |_ a region onii-n plane where c±(u, v) < 

and C2(u, v) > 0. Note that the characteristic speeds are numerated in increasing order, c\ < ci- 



As we can see in FigfJl the splitting contains a special point, an umbilic point, u* = v* = |, 
where characteristic velocities both vanish, c\(u*,v*) = C2(u*,v*) = 0, for any value of Q, as can 
be straighforwardly verified from (j3J) , ([5]) . For Q > Q cr it = — | , the umbilic point is a crossing point 
of the curves c\{u,v) = and C2{u,v) = 0, see FigJJJa). The respective current-density surfaces 
jk(u,v,Q) have a regular convex topology. For Q < Q cr u, current-density surfaces jk(u,v,Q) 
develop a saddle point, and the umbilic point becomes an isolated point, see Figf2]^b). 

Characteristic speeds determine stability of large scale excitations in our system, described on 
the macroscopic scale by a system of conservation laws for coarse- grained densities u(x,t),v(x,t) 

d t u + d x j 1 (u,v) = (7) 
d t v + d x j 2 {u,v) = 0, 



ji and j2 being steady particle currents 



20(|, completemented with boundary conditions 

u(0,t) = u L ; u(l,t) = u R , 
v(0,t) = v L ; v(l,t) = v R . 




Figure 2: Splitting of the physical region into domains G++, G , G |_ according to signs of characteristic 

velocities, for different Q. Boundaries between the domains, on which one characteristic velocity vanishes 
Ci(u,v) = 0, are marked by lines. Point u = v = 1/2 is an umbilic point where c\ = C2 = for any value 
of Q. ( Panel (a): Q > Q C rit- Thick and thin lines stand respectively for Q = —0.5, —0.75. Umbilic pont 
is situated at the crossing of two curves. Panel (b): Q < Q cr it- Lines, in order of increasing thickness, 
correspond to Q = —0.76, —0.8, —0.9, —0.99. Umbilic point, marked by Goo, is an isolated point. 



A commonly made assumption about the flux functions ji,j2, called strict hyperbolicity, reads: 
the characteristic speeds are different c\(u, v) ^ C2{u,v) for all u,v. Strictly hyperbolic systems 



have only two types of fundamental solutions: shocks and rarefaction waves |14fl . Presence of an 
umbilic point ruins strict hyperbolicity and results in appearance of novel excitations listed below. 

An isolated umbilic point in our system (for Q < —3/4) was related to an existence of a large- 
scale excitation, called U-shock: it is microscopically sharp like a shock but according to usual 
stability criteria it should be unstable 



rarefaction waves 



111 ]. The U-shock interface, see Fig|3Ja), is connecting two 



211 ] . If the umbilic point is not isolated (for Q > —3/4 ), the U-shock reduces to 



a bulk homogeneous state with c\ = C2 = 0, meaning that it has densities u* = v* = 1/2, matching 
the umbilic point. We call this state an umbilic state. Both U-shock and umbilic state profiles have 
a property of being left-right symmetric, while the boundary conditions, generically, are not. In 
the next section we describe the domain on the phase diagram, occupied by umbilic point- related 
phases, and the respective phase transitions. 



IV. BOUNDARY-DRIVEN PHASE TRANSITIONS CONTROLLED BY UMBILIC 

POINT 

The splitting of the physical region, done in the previous section, allows to study a phase diagram 
of an open system, and transitions between those, in a systematic manner. A transition from one 
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Figure 3: Average density profiles in phases controlled by umbilic point in a bidirectional model FiglTJ for 
Q = —0.9 < Qcrit{ Panel (a)) and Q = —0.5 > Q C rit( Panel (b)). Symbols are data points obtained from 
Monte Carlo simulations while lines are given by numerical solution of hydrodynamic equations (|B1[) . Values 
of the characteristic velocities are indicated on top of density profiles. Parameters: N = 500, (ul,Vl) = 
(0.9,0.3), (ur,vr) = (0.1,0.9). Note that bulk profiles are left-right symmetric while the boundary rates 
are not. Broken lines on Panel (a) indicate theoretically predicted bulk values, the U-shock amplitude being 
equal to 2A where A is given by ©. Both umbilic point-controlled phases are stable, if (ul,vl) € G__ 
and (u_r,w_r) € G++, see also Tables HHTl 



phase to another in an open system with fixed bulk rates happens by a gradual adiabatic change 
of the boundary rates, which amounts to a respective adiabatic change of boundary densities 
uliVl,ur,vr. Such a gradual change defines a path, or trajectory in the 4-dimensional space 
r(s) = {ul(s),vl(s),ur(s),vr(s)}, parametrized by a running variable s. Choosing the average 
particle densities u, v in the steady state as an order parameter, we study a correspondence T(s) — > 
(u(s),v(s)), which shows singular behaviour at the critical points s cr it along a path. Across a critical 



point a transition between two neighbouring phases takes place, which, see [10| is always governed 
by a large scale excitation- either a shock or a rarefaction wave, depending on how boundary 
changes are performed. To specify, denote the number of positive characteristic velocities in the 
left (right) boundary reservoir as #l (#_r), which can then take integer values #_r,#l = 0, 1,2. 
A transition from one phase to another is governed by a rarefaction wave if #l(s) < #r(s) for all 
s. A straightforward example of a path T(s) satisfying the requirement #l(s) < #_r(s), is a Ph2 



path described by steps I, II, III below, see also first two columns of Tab 
system a mapping T(s) — > (u(s),v(s)) along a Ph2 path is continuous 

An isolated umbilic point, present in our model for large interlane interactions Q < —3/4, 



.ejp In a strictly hyperbolic 



stabilizes a sharp interface connecting two rarefaction waves, a so-called a U-shock 



111 ]. Below we 



demonstrate that the U-shock makes possible a discontinuous change T(s) — > (u(s),v(s)) along a 
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Ph2 path in a system with an isolated umbilic point. 

In order to clarify the influence of a U-shock on the phase diagram we consider a Ph2 path 



2 2[ | from a steady state G++ to a steady state G , in presence and in absence of an isolated 

umbilic point. Let us denote by < s < 1 a variable parametrizing the adiabatic Ph2 path, and 
by ul(s),vl(s),ur(s),ur(s) the respective densities of boundary reservoirs. As in [lOj], we shall 
vary the boundary densities along the Ph2 path in the following way: 

I. Initial point s = and final point s = 1 corresponds to fully-matching left and right bound- 
ary reservoirs. u L (0) = u R {0) = u initial ,v L (0) = v R (0) = v initial , where Cl (u initial jV initial ) > 
Q^^nitial ynitial^ > _ Analogously for u L (l) = u R (l) = u final , v L (0) = v R (0) = v final , where 
Cl {u final ,v final ) < 0,c 2 (u final ,v final ) < 0. 

II. For < s < ~, left boundary densities ul(s), vl(s) are changing smoothly from u initial ; v imtiai 
to u f mal jV f mal at s = 1/2. The right boundary densities remain the same, u R (s < 1/2) = 

III. For 2 < s < 1, the right boundary densities u R (s),v R (s) are changing smoothly from 
u imtiai ^initial g _ u f mal ^ y f mal at s _ ^ -p^g } e fl- boundary densities remain the same, 

u R (s > 1/2) = uf inal ,v R (s > 1/2) = v final . 

As argued in all along such a path the steady state is controlled by rarefaction waves. 
In absence of an umbilic point, rarefaction waves do not have any discontinuities [l^.fliil]. and 
therefore the stationary densities along the path u(s),v(s) are expected to change continuously 
with s. 

In the following we demonstrate that an isolated umbilic point provokes two discontinuous 
jumps of u(s),v(s), and locate the critical point. The jumps are due to a motion of the U-shock 
between the boundaries which become biased at the transition point. 

We choose the initial and the final state to be ( u initial jV initial ) = (0.1,0.9) and (uf inal , v f inal ) = 
(0.9, 0.3), and change the boundary densities for intermediate s by linear interpolation, i.e. ul(s) = 
u initial + 2 s (u final - u initial ), v L (s) = v initial + 2s(v^ inal - v initial ) for < s < 1/2, and similarly 
for s > 1/2. After that we perform an adiabatic Ph2 path for a system with an isolated umbilic 
point Q = —0.9 < Qcrit, and non-isolated umbilic point Q = —0.5 < Q C rit- For both values of Q, 

^initial ^initial} e q ++ and ( u final ^finalj £ Wg present the resU ltS ill Figs. SEE 

Case 1. Isolated umbilic point. An isolated umbilic point appears above the critical 
interaction amplitude, Q < — |. In this case, Fig. IU(a) , we see a discontinuous change in average 
stationary densities u(s) and v(s) at two points sl rit ,s1 rit . As already stressed, this discontinuous 
phase transition is impossible in a strictly hyperbolic system. A closer examination reveals the 
mechanism of the new transition. Firstly, we notice that the steady state of a system is a U-shock 
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Figure 4: Panel (a): Location of stationary bulk densities u(s), v(s) along the Ph2 path on a u — v plane, in 
case of isolated umbilic point. Points A and B mark initial and final point of the path, respectively. Curves 
marked Go+ and G_o denote location of points with ci(u, v) = and 02(11, v ) = 0, respectively. Broken line 
mark a discontinuous transition to the U-shock phase. Spreading of the points in the middle is due to finite 
size-effects. Panel (b): Stationary bulk densities u(s), v(s) of the right movers (small circles) and of the left 
movers (triangles) along the Ph2 path. Broken lines denote the dependence of (ul{s),vl(s)) for s < 1/2 and 
( u r(s): v r{ s )) f° r s > 1/2- Thin lines denote Ck(u,v). Critical points s cr mark discontinuous transitions to 
the U-shock phase. Most data arc obtained bv numerical integration of the mcanficld equations. Large circles 
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Figure 5: Panel (a): Location of stationary bulk densities u(s), v(s) along the Ph2 path onaa - v plane, 
for Q = —0.5. Points A and B mark initial and final point of the path, respectively. Curves marked Go+ 
and G-o denote location of points with ci(u,v) = and C2(u,v) = 0, respectively. Panel (b): Stationary 
bulk densities u(s),v(s) of the right movers (circles) and of the left movers (triangles) along the Ph2 path. 
Broken lines denote the dependence of («l(s), vl(s)) for s < 1/2 and (ur(s), Vr(s)) for s > 1/2. Thin lines 
denote cu{u,v). Critical points s l crit and s\ rit mark transitions to the umbilic phase. Results are obtained 
by numerical integration of the meanfield equations. Large circles on Panel(b) mark steady state densities 
obtained from Monte-Carlo simulations. Parameters: O — —0.5, N = 200. 
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if 

(u L ,v L )eG ++ and (ur, vr) G G — , (8) 

which corresponds to a segment s G [slriti s crit] °f * ne adiabatic Ph2 trajectory. The U-shock is 
an interface between the two rarefaction waves. The average bulk densities to the left and to the 
right from the interface depend only on Q and are given by u, v = \ ± A, where 

A = iv / 3Q =I +4. (9) 

The above values of u,v correspond to maxima of the current density relation j(u,v,Q), which 
develops a saddle point for Q < —3/4, see Figj6j Steady-state currents of both species have the 
same amplitude, f{ = = 1/(8 |Q|). 

As long as the conditions ([8]) are satisfied, the bulk densities, and consequently the steady 
currents, do not depend on boundary densities. The boundary densities affect only the respective 
boundary layers, interpolating between the bulk and the boundary, on a microscopic scale. The 
interface position fluctuates with time around the center of the lattice. 

Now consider a left vicinity of a critical point s = s\ rit — e, where < e < 1. Once e > 0, the 
U-shock is biased to the left, and gets pinned to the left boundary. The resulting steady state is 
homogeneous with the densities u = h — A,v = h + A. With e increasing, (u(s),v(s)) follows the 
curve G + q. 

Analogously, at the other side of the U-shock phase, in the right vicinity of another critical point 
s = s 2 crit + s, the U-shock is biased to the right. This results in a homogeneous steady state with 
the densities u = ^ + A,v = ^ — A, where A is given by ([9]). As s increases, the point (u(s),v(s)) 
follows the curve Go-. Depinning from the curve Go- occurs when both (ul,vl) G G ++ and 
(ur,v r ) G G ++ . 

Steady states and corresponding reservoir densities are listed in order of their temporal appear- 
ance along the adiabatic trajectory on Table HI As one can see, steady states are straightforwardly 
connected to the phase space splitting in subregions G Q( g with different signs of characteristic 
velocities. 

It is instructive to compare the discontinuous U-shock-governed phase transition described above 



to a usual discontinuous phase transition in driven systems governed by a standard Lax shock 13J |. 
In the latter, the shock changes the sign of a bias at the critical point. Exactly at the critical point, 
the shock is unbiased, and, due to fluctuations, performs a random walk between the boundaries. 
On the contrary, in a U-shock -controlled phase transition, the U-shock stays unbiased in the 
whole segment [sl rit , s^, rit ] , defined by (jSJ). For s within this segment, the U-shock position is not 
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(ul,vl) 


{ur,vr) 


Steady state 


Steady state 


belong to 


belong to 


belongs to 


densities 


G ++ 


G ++ 


G++ 


U = U L , V = V L 


G-+ 


G ++ 


Gq+ 


pinned to Go+ 


G 


G++ 


G0+/G-0 


U-shock u = 5 ± A, u = i =F A 


G 


G-+ 


G_ 


pinned to G_o 


G 


G__ 


G__ 


U = U R , V = V R 



Table I: Steady states densities u(s), v(s) and corresponding particle densities in boundary reservoirs along 
a Ph2 path. A is given by ©. 



{ul,v l ) 


(ur,vr) 


Steady state 


Steady state 


belong to 


belong to 


belongs to 


densities 


G__ 


G++ 


Goo 


umbilic point rarefaction u = v = 5 



Table II: Open bidirectional system with a non-isolated umbilic point. Steady states and corresponding 
boundary reservoir densities along the middle part of a Ph2 path. For remaining path segments, see Table Q] 

performing a random walk between the boundaries, but fluctuates around the middle point as if it 
was in a potential well. 

Case 2. Non-isolated umbilic point. Now, let us push the interaction Q towards the 
critical point Q cr it = —3/4. The U-shock amplitude 2 A becomes smaller and smaller until it 
disappears at the critical point A(Q cr it) = 0. The whole U-shock-governed phase above the 
critical point Q > Q c nt reduces to a homogeneous phase with densities matching the umbilic point 
c\(u, v) = C2(u,v) = 0, see Table HU In this way, the umbilic point defines a robust phase on the 
phase diagram, which we shall call an umbilic phase or U-phase. The U-phase with c\ = C2 = 

appears to be stable whenever (ul,vl) € G and (ur,vji) G G ++ , see Table HU in the very 

same domain where a U-shock is stable ([8]). A pinning-depinning transition from the U-phase is 
continuous, in contrary to a transition from/to U-shock -governed phase for Q < Q cr u discussed 
earlier, see Figj5j 

Several comments are in order at this point. Firstly, note that the steady state particle currents 
and bulk density profiles for an umbilic phase and for a U-shock phase are left-right symmetric (so 
as the bulk hopping rates), in spite of boundary conditions being explicitly not left-right symmetric. 
Away from the umbilic (or U-shock) phase, a steady state is sensitive to the boundaries and in 
general is not left-right symmetric. In this way, when entering the respective U-phases, the steady 
state becomes insensitive to boundaries and regains its bulk symmetry. In absence of umbilic points 
(i.e. in strictly-hyperbolic systems), steady states are generically boundary-rates dependent. 
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Figure 6: Cuts of the surface j\{u,v) along the lines u = v (the convex curve) and u = 1 — v (the curve 
with double maximum), for Q < Q cr it- Parameters: Q = —0.9. Points u = 1/2 ± A correspond to current 
maxima, A is given by ([9]). Above the critical point Q > Q cr it both cuts are convex. 

Secondly, we observe that the total current of both species j\ + \ along the Ph2 path attains 
its maximal value just for an umbilic phase and of a U-shock phase. This fact exemplifies a validity 
of the maximal current principle for bidirectional systems. Indeed, in systems with one driven 
particle species and open boundaries the maximal current principle asserts in particular that for 
maximal feeding regime (when both entrance and exit rates are maximal), the stationary current 
is maximized with respect to an average particle density p, j steady = max pe[o,p max ] j(p)- I n our 
model with two species, such a maximal feeding regime is realized when (ul,vl) = (1,0) and 
(ur, vr) = (1,0), which, consulting the Fig|2]and ([8]), corresponds to the U-shock or Umbilic phase 
domain, depending on the value of Q. The respective maximal value of the total current is j\ + \j2 \ = 
2ji = (4 IQir 1 for U-shock Q < Q crit and j x + |j 2 | = 2j x (u = \,v = \) = ^Q + l/{VQ + ^ + 1) 
for the umbilic phase Q > Qcrit, as can be straighforwardly derived from analytic expressions for 
the currents. 

V. CONCLUSIONS 

We conclude that a presence of an umbilic point in the current-density function, isolated or not, 
gives rise to new types of boundary driven phase transitions. If the umbilic point is an isolated 
one, one has a discontinuous transition from Go+ steady state to the U-shock state Go+ /G_o and 
another discontinuous transition from the U-shock state to G_o state, both transitions governed 
by a biased motion of a U-shock. In case of non-isolated umbilic point, two continuous (pinning- 
depinning) transitions Go+ — > Goo an d Goo ~~ ^ G-o take place. All these transitions are observable 
along any Ph2 path, leading from a state with positive characteristics to a steady state with 
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negative characteristics via adiabatically changing boundary conditions. There is no hysteresis 
of any kind, so by inversion of the path the sequence of transitions is inverted. The U-shock 
phase and umbilic phases are robust, and exemplify the validity of a maximal current principle for 
bidirectional particle models. We identified the domain on the phase diagram occupied by umbilic 
point related phases (J8]). Within this domain the system regains its bulk left-right symmetry in 
spite of the boundary conditions being explicitly not left-right symmetric. 

How robust are our results? First of all, the umbilic point with c\ = C2 = is a general feature of 
a models with left-right symmetry of hopping rates U(], of which we considered a special example 



with an exactly known steady state on a ring. We expect qualitatively similar results for general 
particle systems whether "solvable" or not. In an open system it is only necessary to maintain a 
stationary maximal flow regime to find out whether the umbilic point is isolated or not. Appearance 
of a discontinuity in a bulk density profile in the maximal flow regime would indicate that an umbilic 
point has become isolated. In such a way from a simple single macroscopic observation of a system 
one can judge on intrinsic differential properties of its current-density function. 

Bidirectional models are being widely studied in the literature, in particular, in connection with 
the intriguing phenomenon of spontaneous symmetry breaking (SSB) [23|]-32[|. In most studies 
however the current-density function is a convex surface. It would be interesting to study SSB 
in presence of an isolated umbilic point. It is interesting to note that bidirectional models have 
also an integrability aspect 3], however up to now no example of an integrable system with open 
nontrivial boundary conditions have been presented. 
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Appendix A: Boundary rates 

Boundary rates for injection and extraction of the particles are chosen so that constant densities 
of particles are kept on the left and on the right boundary. Steady state of our model with the 
rates ([I]) for a system on a ring has a remarkable property: for a configuration C the stationary 
probability is given by a product measure 

Pc = Z- 1 J]e~ imkmk , (Al) 

k 
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where n k ,m,k = 0,1 are particle occupation number on site k on chains 1 and 2, and Z is a 
normalization. We see that neighbouring pairs of adjacent sites are uncorrelated. This fact allows 
to express all steady state equal-time expectations in terms of probabilities of a single pair of 
adjacent sites £l nkmk (u,v) where u, v is an average density of right and left movers. We have 
Qqo + Qu + Q,qi + Oio = 1, J^io + On = u, Qqi + = v -> an d v) is given by (|6|). Procedure, 



completely analogous to that in 33J] results in the following definition of the boundary rates: a 
right-moving particle is injected to a site 1 with rate ul (with rate ul + QQu(ul,vl) ) if an 
adjacent site is empty (filled). A right-moving particle is extracted from site N with rate 1 — ur 
(with rate 1 — ur + Q^oo{ur,vr) ) if an adjacent site is filled (empty). Analogously, left moving 
particles are injected at the right boundary with rate vr (with rate vr + Q£1±i(vr, ur) ) if an 
adjacent site is empty (filled), and extracted at the left boundary with rate 1 — vl (with rate 
1 — vl + QQoo(vl,vl) ) if an adjacent site is filled (empty). In case of matching left and right 
boundaries, ul = ur = u, vl = vr = v, the exact steady state of the system is (1A1|) . independently 
of system size N. 

Appendix B: Meanfield equations 

In our meanfield approximation, we neglect correlations between the adjacent pairs of sites, 
which are also absent in the steady state (|Aip but not between the adjacent sites. The equations 
are obtained by averaging the exact operator equations of motion for occupation number operators 
f^ki^kt which for right-moving particle at site k read 

d{n k ) r > ,„ > 

Q t = \3k-l,k) - \Jk,k+l) 

where jk,k+i = {n k (l-h k +i)) + Q{n k m k (l-h k+ i)(l-m k+1 )) . Denoting {h k (t)) = s k (t), {m k (t)) = 
q k (t), we simplify parts of above expression as (h k (l — %+].)) ~ {n k )(l — n k+ \) = s k {l — q k +\) 
and (n k m k (l - h k+1 )(l - rh k+1 )) w (n k m k )((l - n k+x )(l - rh k+1 )) = flu(s k ,q k )n 00 (s k+1 ,q k+1 ), 
where £ln(u, v) is given by exact microscopic expression ©. Thus, for a homogeneous state 
{n k ) = u i (iTT-k) = v the meanfield expression gives the exact microscopic stationary current @. 
The semiclassical equation of motion becomes 

ds k 

— = +Sfe_i(l - q k ) - QOn(sfc-i,gfc-i)0 o(sfc,g , fc) 

- s k (l - q k+ i) - Qtl u (s k , q k )Q 00 (s k+ i,q k+ i), (Bl) 

for k = 2,3, ...N — 1, complemented with boundary conditions for the boundary sites k = 1 and 
k = N 

si(t) = u L ; s N (t) = u R ; (B2) 
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For left-movers the equations are derived analogously. The complete set of equations of motion can 
be viewed as a natural discretization scheme with which we integrate numerically the hydrodynamic 
equations ([7]). Indeed, by Taylor expansion of (|B1|) . and Euler rescaling of space and time we obtain 
(f7]). Comparison with the stochastic evolution shows that both steady state profiles and temporal 
evolution is described correctly. 
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